#include "RefSpec.h"
#include <assert.h>


#define LOG2(x) log2(x)//(log(double(x))/log(double(2)))
typedef std::vector<double> Array1d;

void CalRefPSD(const REF_SEGMENT& sect, const Array1d& vf,
	Array1d& vPsd,
	Array1d& vAlarm1, Array1d& vAlarm2,
	Array1d& vAbort1, Array1d& vAbort2);
void CalRefCoh(const REF_SEGMENT& sect, const Array1d& vf,
	Array1d& vCohPhi,
	Array1d& vAlarm1, Array1d& vAlarm2,
	Array1d& vAbort1, Array1d& vAbort2);
void CalRefPhi(const REF_SEGMENT& sect, const Array1d& vf,
	Array1d& vPhi,
	Array1d& vAlarm1, Array1d& vAlarm2,
	Array1d& vAbort1, Array1d& vAbort2);


void CalRefPSD(const REF_SEGMENT& sect, const Array1d& vf,
	Array1d& vPsd,
	Array1d& vAlarm1, Array1d& vAlarm2,
	Array1d& vAbort1, Array1d& vAbort2) 
{
	assert(vf.size() == vPsd.size());

	double df = sect.df;
	auto psd1 = sect.val1;
	auto psd2 = sect.val2;
	auto slope = sect.slope;

	int type = 0;
	if (!slope.empty() && !psd2.empty())// [Nan, m, psd2]
	{
		type = 1;
	}
	else if (!psd1.empty()&& !psd2.empty())// [psd1, Nan, psd2]
	{
		type = 2;
	}
	else if (!psd1.empty() && !slope.empty())// [psd1, m, Nan]
	{
		type = 3;
	}

	double PSD1 = atof(psd1.data());
	double PSD2 = atof(psd2.data());
	double m = atof(slope.data());

	double f1 = sect.f1;
	double f2 = sect.f2;

	int n1 = int(f1 / df);
	int n2 = int(f2 / df);

	auto& R = vPsd;
	auto& f = vf;

	switch (type)
	{
	case 1: // [Nan, m, psd2]
		for (int i = n1; i<=n2; i++)
		{
			R[i] = PSD2*pow(10., 0.1*m*LOG2(f[i] / f2));
		}
		break;
	case 2:// [psd1, Nan, psd2]
		m = 10 * log10(PSD2 / PSD1) / LOG2(f2 / f1);
		for (int i = n1; i <= n2; i++)
		{
			R[i] = PSD2*pow(10., 0.1*m*LOG2(f[i]) / f2);
		}
		break;
	case 3:// [psd1, m, Nan]
		for (int i = n1; i <= n2; i++)
		{
			R[i] = PSD1*pow(10., (0.1*m*LOG2(f[i] / f1)));
		}
		break;
	}

	double t0 = 0, t1 = 0;
	double db = 0;
	for (int i = n1; i <= n2; i++)
	{
		t0 = R[i];

		db = atof(sect.alarm1.data());
		t1 = (log10(t0) * 10 + db) / 10;
		vAlarm1[i] = pow(10., t1);

		db = atof(sect.alarm2.data());
		t1 = (log10(t0) * 10 + db) / 10;
		vAlarm2[i] = pow(10., t1);

		db = atof(sect.abort1.data()); 
		t1 = (log10(t0) * 10 + db) / 10;
		vAbort1[i] = pow(10., t1);

		db = atof(sect.abort2.data());
		t1 = (log10(t0) * 10 + db) / 10;
		vAbort2[i] = pow(10., t1);
	}
}

void CalRefCoh(const REF_SEGMENT & sect, const Array1d & vf, Array1d & vCohPhi, Array1d & vAlarm1, Array1d & vAlarm2, Array1d & vAbort1, Array1d & vAbort2)
{
	assert(vf.size() == vCohPhi.size());

	double df = sect.df;
	auto val1 = sect.val1;
	auto val2 = sect.val2;
	auto slope = sect.slope;

	int type = 0;
	if (!slope.empty() && !val2.empty())// [Nan, m, psd2]
	{
		type = 1;
	}
	else if (!val1.empty() && !val2.empty())// [psd1, Nan, psd2]
	{
		type = 2;
	}
	else if (!val1.empty() && !slope.empty())// [psd1, m, Nan]
	{
		type = 3;
	}

	double dVAL1 = atof(val1.data());
	double dVAL2 = atof(val2.data());
	double m = atof(slope.data());

	double f1 = sect.f1;
	double f2 = sect.f2;

	int n1 = int(f1 / df);
	int n2 = int(f2 / df);

	auto& R = vCohPhi;
	auto& f = vf;

	switch (type)
	{
	case 1: // [Nan, m, psd2]
		for (int i = n1; i <= n2; i++)
		{
			R[i] = dVAL2 - (f2 - f[i])*m;
		}
		break;
	case 2:// [psd1, Nan, psd2]
		m = (dVAL2 - dVAL1) / (f2 - f1);
		for (int i = n1; i <= n2; i++)
		{
			R[i] = dVAL2 - (f2 - f[i])*m;
		}
		break;
	case 3:// [psd1, m, Nan]
		for (int i = n1; i <= n2; i++)
		{
			R[i] = dVAL1 + (f[i] - f1)*m;
		}
		break;
	}

	double t0 = 0, t1 = 0, db = 0;
	for (int i = n1; i <= n2; i++)
	{
		t0 = R[i];

		db = atof(sect.alarm1.data()) / 100;
		vAlarm1[i] = t0*(1+db);

		db = atof(sect.alarm2.data()) / 100;
		vAlarm2[i] = t0*(1 + db);

		db = atof(sect.abort1.data()) / 100;
		vAbort1[i] = t0*(1 + db);

		db = atof(sect.abort2.data()) / 100;
		vAbort2[i] = t0*(1 + db);

	}
}

void CalRefPhi(const REF_SEGMENT & sect, const Array1d & vf, Array1d & vPhi, Array1d & vAlarm1, Array1d & vAlarm2, Array1d & vAbort1, Array1d & vAbort2)
{
	assert(vf.size() == vPhi.size());

	double df = sect.df;
	auto val1 = sect.val1;
	auto val2 = sect.val2;
	auto slope = sect.slope;

	int type = 0;
	if (!slope.empty() && !val2.empty())// [Nan, m, psd2]
	{
		type = 1;
	}
	else if (!val1.empty() && !val2.empty())// [psd1, Nan, psd2]
	{
		type = 2;
	}
	else if (!val1.empty() && !slope.empty())// [psd1, m, Nan]
	{
		type = 3;
	}

	double dVAL1 = atof(val1.data());
	double dVAL2 = atof(val2.data());
	double m = atof(slope.data());

	double f1 = sect.f1;
	double f2 = sect.f2;

	int n1 = int(f1 / df);
	int n2 = int(f2 / df);

	auto& R = vPhi;
	auto& f = vf;

	switch (type)
	{
	case 1: // [Nan, m, psd2]
		for (int i = n1; i <= n2; i++)
		{
			R[i] = dVAL2 - (f2 - f[i])*m;
		}
		break;
	case 2:// [psd1, Nan, psd2]
		m = (dVAL2 - dVAL1) / (f2 - f1);
		for (int i = n1; i <= n2; i++)
		{
			R[i] = dVAL2 - (f2 - f[i])*m;
		}
		break;
	case 3:// [psd1, m, Nan]
		for (int i = n1; i <= n2; i++)
		{
			R[i] = dVAL1 + (f[i] - f1)*m;
		}
		break;
	}

	double t0 = 0, t1 = 0, db = 0;
	for (int i = n1; i <= n2; i++)
	{
		t0 = R[i];

		db = atof(sect.alarm1.data());
		vAlarm1[i] = t0 + db;

		db = atof(sect.alarm2.data());
		vAlarm2[i] = t0 + db;

		db = atof(sect.abort1.data());
		vAbort1[i] = t0 + db;

		db = atof(sect.abort2.data());
		vAbort2[i] = t0 + db;
	}
}



void REF_DATA::ResetVector()
{
	vf.resize(nFFT / 2 + 1, 0);
	vRef.resize(nFFT / 2 + 1, 0);
	vAlarm1.resize(nFFT / 2 + 1, 0);
	vAlarm2.resize(nFFT / 2 + 1, 0);
	vAbort1.resize(nFFT / 2 + 1, 0);
	vAbort2.resize(nFFT / 2 + 1, 0);

	//vf.clear();
	//vRef.clear();
	//vAlarm1.clear();
	//vAlarm2.clear();
	//vAbort1.clear();
	//vAbort2.clear();
}

void REF_DATA::CalOverallCoh()
{
	ResetVector();
	if (listSegment.empty())
	{
		return;
	}
	double df = listSegment.front().df;
	for (int i = 0; i < vf.size(); i++)
	{
		vf[i] = i*df;
	}
	for (auto it = listSegment.begin(); it != listSegment.end(); it++)
	{
		auto& segment = *it;
		::CalRefCoh(segment, vf, vRef, vAlarm1, vAlarm2, vAbort1, vAbort2);
	}
}

void REF_DATA::CalOverallPhi()
{
	ResetVector();
	if (listSegment.empty())
	{
		return;
	}
	double df = listSegment.front().df;
	for (int i = 0; i < vf.size(); i++)
	{
		vf[i] = i*df;
	}
	for (auto it = listSegment.begin(); it != listSegment.end(); it++)
	{
		auto& segment = *it;
		::CalRefPhi(segment, vf, vRef, vAlarm1, vAlarm2, vAbort1, vAbort2);
	}
}



void REF_PSD::CalOverallPSD()
{
	ResetVector();
	if (listSegment.empty())
	{
		return;
	}
	double df = listSegment.front().df;
	for (int i = 0; i < vf.size(); i++)
	{
		vf[i] = i*df;
	}
	for (auto it = listSegment.begin(); it != listSegment.end(); it++)
	{
		auto& segment = *it;
		::CalRefPSD(segment, vf, vRef, vAlarm1, vAlarm2, vAbort1, vAbort2);
	}

	// CalRms
	int n1 = int(listSegment.front().f1 / df);
	int n2 = int(listSegment.back().f2 / df);
	dRMS = 0;
	for (int i = n1; i <= n2; i++)
	{
		dRMS += ::sqrt(vRef[i] * df);
	}
	// Max Acc
	this->dMaxAcc = 3 * dRMS;

}

void REF_PSD::RescaleRMS(double newRMS)
{
	double r = ::pow(newRMS/ dRMS, 2);
	double newPSDVal = 0;
	for (auto it = listSegment.begin(); it != listSegment.end(); it++)
	{
		char buf[20];

		memset(buf, 0, sizeof(buf));
		if (!it->val1.empty())
		{
			newPSDVal = ::atof(it->val1.data())*r;
			::sprintf_s(buf, sizeof(buf), "%g", newPSDVal);
		}
		it->val1 = buf;

		memset(buf, 0, sizeof(buf));
		if (!it->val2.empty())
		{
			newPSDVal = ::atof(it->val2.data())*r;
			::sprintf_s(buf, sizeof(buf), "%g", newPSDVal);
		}
		it->val2 = buf;
	}
	CalOverallPSD();
}

void REF_DATA::WriteFile(const std::string & szFile)
{
	std::ofstream outfile;
	outfile.open(szFile.data());
	outfile << *this << std::endl;
}


std::ostream & operator<<(std::ostream & os, const REF_DATA & Ref)
{
	os << Ref.szName << std::endl;
	os << Ref.nFFT<< std::endl;
	os << Ref.listSegment.size() << std::endl;
	for (auto it = Ref.listSegment.begin(); it != Ref.listSegment.end(); it++)
	{
		os << it->f1 << std::endl;
		os << it->f2 << std::endl;
		os << it->df << std::endl;
		os << it->val1 << std::endl;
		os << it->slope << std::endl;
		os << it->val2 << std::endl;
		os << it->alarm1 << std::endl;
		os << it->alarm2 << std::endl;
		os << it->abort1 << std::endl;
		os << it->abort2 << std::endl;
	}
	return os;
}

std::istream & operator >> (std::istream & is, REF_DATA & Ref)
{
	char buf[100] = { '\0' };
	is.getline(buf, sizeof(buf));
	Ref.szName = buf;
	is.getline(buf, sizeof(buf));
	Ref.nFFT = ::atoi(buf);

	is.getline(buf, sizeof(buf));
	int n = ::atoi(buf);
	Ref.listSegment.resize(n);
	for (auto it = Ref.listSegment.begin(); it != Ref.listSegment.end(); it++)
	{
		REF_SEGMENT segment;
		is.getline(buf, sizeof(buf));
		segment.f1 = atof(buf);
		is.getline(buf, sizeof(buf));
		segment.f2 = atof(buf);
		is.getline(buf, sizeof(buf));
		segment.df = atof(buf);
		is.getline(buf, sizeof(buf));
		segment.val1 = buf;
		is.getline(buf, sizeof(buf));
		segment.slope = buf;
		is.getline(buf, sizeof(buf));
		segment.val2 = buf;
		is.getline(buf, sizeof(buf));
		segment.alarm1 = buf;
		is.getline(buf, sizeof(buf));
		segment.alarm2 = buf;
		is.getline(buf, sizeof(buf));
		segment.abort1 = buf;
		is.getline(buf, sizeof(buf));
		segment.abort2 = buf;
		*it = segment;
	}
	is.getline(buf, sizeof(buf));


	return is;
}
